% clear all
% %clc
%
%load('data_exp.mat')

N_temp = 500;

% combos_all_unique = combos_all;
combos_all_unique = double(combos_all_unique);

% [data,data_trans] = generate_data_rand(N_temp,J,mu_x,2,mu_z,2,alpha,beta,fraction);

ratio_temp_second = zeros(N_temp,1);
ratio_temp_first = zeros(N_temp,1);
d2s_dzdz = zeros(N_temp,1);
d2s_dzdx = zeros(N_temp,1);
ds_dz = zeros(N_temp,1);
ds_dx = zeros(N_temp,1);
ratio_trimmed_mean_temp_second = NaN(2,1);
ratio_trimmed_mean_temp_first = NaN(2,1);
trimmed_mean_of_ratios_temp = NaN(2,1);

for j=1:2
    
    parfor i=1:N_temp

        % second derivatives
        
        [~,temp] = max(data.z(i,:));
        
        temp_diff = setdiff(1:J,temp);
        
        if j==2
            temp_diff = temp_diff([2 1]);
        end
        
        ind_good2 = temp_diff(1);

        arg_x_1 = normcdf(data.x(i,temp));
        arg_x_other = normcdf(data.x(i,temp_diff));

        arg_z_1 = normcdf(data.z(i,temp));
        arg_z_other = normcdf(data.z(i,temp_diff));

        argument = [arg_x_1 arg_x_other arg_z_1 arg_z_other];

        deriv_z_1 = normpdf(data.z(i,temp));
        deriv_z_2 = normpdf(data.z(i,ind_good2));
        deriv_x_2 = normpdf(data.x(i,ind_good2));

        d2s_dzdz(i) = fun_d2sigma_new2(theta_NPD,argument,m_regr,combos_all_unique,J,J+1,J+2)*deriv_z_1*deriv_z_2;

        d2s_dzdx(i) = fun_d2sigma_new2(theta_NPD,argument,m_regr,combos_all_unique,J,J+1,2)*deriv_z_1*deriv_x_2;

        ratio_temp_second(i) = d2s_dzdz(i)/d2s_dzdx(i);


        % first derivatives

        argument = normcdf([data.x(i,:) data.z(i,:)]);

        deriv_x_1 = normpdf(data.x(i,1));
        deriv_z_1 = normpdf(data.z(i,1));

        ds_dz(i) = fun_dsigma(theta_NPD,argument,m_regr,combos_all_unique,J+1)*deriv_z_1;

        ds_dx(i) = fun_dsigma(theta_NPD,argument,m_regr,combos_all_unique,1)*deriv_x_1;

        ratio_temp_first(i) = ds_dz(i)/ds_dx(i);


    end
    
    [num_temp_ord,ind_num] = sort(d2s_dzdz);
    [den_temp_ord,ind_den] = sort(d2s_dzdx);
    ratio_trimmed_mean_temp_second(j) = mean(num_temp_ord(N_temp*0.25:N_temp*0.75))/mean(den_temp_ord(N_temp*0.25:N_temp*0.75));

    [num_temp_ord,ind_num] = sort(ds_dz);
    [den_temp_ord,ind_den] = sort(ds_dx);
    ratio_trimmed_mean_temp_first(j) = mean(num_temp_ord(N_temp*0.25:N_temp*0.75))/mean(den_temp_ord(N_temp*0.25:N_temp*0.75));

    [ratio_temp_ord,ind] = sort(ratio_temp_second);
    trimmed_mean_of_ratios_temp(j) = mean(ratio_temp_ord(N_temp*0.25:N_temp*0.75));
    
    %i
end

ratio_trimmed_mean_second = mean(ratio_trimmed_mean_temp_second)
ratio_trimmed_mean_first = mean(ratio_trimmed_mean_temp_first)
trimmed_mean_of_ratios = mean(trimmed_mean_of_ratios_temp)

% ratio = median(ratio_temp);
% 
% ratio_means = mean(d2s_dzdz)/mean(d2s_dzdx);
% ratio_medians = median(d2s_dzdz)/median(d2s_dzdx);
% 
% 
% [ratio_temp_ord,ind] = sort(ratio_temp);
% ratio_trimmed = mean(ratio_temp_ord(N_temp*0.10:N_temp*0.90));
% 
% regr_coeff = regress(d2s_dzdz,d2s_dzdx); % regr_coeff = regress(d2s_dzdz,[ones(length(d2s_dzdz),1) d2s_dzdx])
% 
% slope = regr_coeff(1); %regr_coeff(2);

% mean(slope)
% std(slope)
%
%
% mean(ratio_means)
% std(ratio_means)

%
% figure(1)
% scatter(d2s_dzdz,d2s_dzdx)
% axis([-1 0 -1 0])
% refline(1,0)
%
%
% [b,bint] = regress(d2s_dzdz, d2s_dzdx)